Repository navigation
Conversation
…ble weights in twas_weights_conversion
gwas_sumstats_construct.R writes no GwasSumStats file for a block without GWAS variants, so twas.R failed on readRDS and the whole region errored. Such genes now get an empty result with a noResult attribute.
twas_manifest.R gains --unit gene (one unit per gene spanning its weight window, so GWAS and LD are read per gene) next to the existing region grouping, which stays the script default. An optional CV prefilter (--min-cv-rsq / --max-cv-pval) skips genes where no context/method passes, so their GWAS input is never built; weights without per-row CV results are kept for pecotmr to decide. twas_ctwas.ipynb uses the gene unit and the prefilter by default (--twas_unit, --twas_prefilter).
- twas_ctwas.ipynb keeps region as the default TWAS unit, so existing commands give the same outputs; the gene unit is opt-in with --twas_unit gene. - CV prefilter: off by default in the notebook (--twas_prefilter turns it on for twas runs), since the ctwas steps read the same manifest and cTWAS uses every gene with weights. It now reads the per-method CV results of the pecotmr 0.9.5+ TwasWeights layout (it previously found no metrics there and kept every gene), ignores methods without a CV r2 such as ensemble, uses the twas step's --rsq-option and --rsq-pval-option, keeps genes with no CV metrics, no cvResult column or an unreadable weight file (reported on stderr), treats --min-cv-rsq <= 0 as off, and loads pecotmr before reading the weights. - Gene unit: make region_id file-name safe and stop when the xQTL meta repeats a gene, since both would collide on output paths. - twas_weights_conversion.R: skip all-zero/NA weights before resolving variant IDs, and prefer preset_variants_result IDs when their length fits. - Tests: the notebook twas test pins the region unit; new tests cover the gene unit (same result as the region unit on the toy gene), gene-unit manifest rows, repeated genes, and the prefilter keep/drop rules.
- twas_manifest.R: --unit gene stops early when a gene has no start/end in the xQTL meta, checks for repeated genes before making IDs file-name safe, and stops when two genes would share a file name. The CV prefilter skips cvResult entries that are not lists instead of reporting the file as unreadable. - twas.R: a missing GwasSumStats file still means an empty window, but a missing folder now stops as a wrong path; the noResult text no longer assumes an LD block. - twas_ctwas.ipynb: the twas step removes an old GwasSumStats file before rebuilding it, makes the gene part of output names file-name safe, and the Input and Command Interface sections list --twas_unit, --twas_prefilter, --rsq_option and --rsq_pval_option. - twas_weights_conversion.R: report genes without usable weights with message(), since the warning was lost on quit(). - Tests: prefilter p-value cutoff (keep and drop, with the pval fall-back and ensemble skip), NA gene window, file-name safe IDs and clashes, missing GwasSumStats file and folder, and conversion of weights without names, with unmatched IDs and with all-zero weights.
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Five commits needed to run genome-wide TWAS with the converted FunGen-xQTL weights on pecotmr 0.9.6. They have been used for the MSBB eQTL (17,297 genes) and KNIGHT pQTL runs on the 2026 AD GWAS.
twas_weights_conversion.R: handle legacy weights stored for the QC'd variant subset only, and genes with no usable weights, instead of failing the conversion.twas.R,twas_manifest.R:gwas_sumstats_construct.Rwrites no GwasSumStats file for an LD block without GWAS variants, sotwas.Rfailed onreadRDSand the whole region errored. Genes there now get an empty result with anoResultattribute.twas_manifest.R,twas_ctwas.ipynb: new--unit gene(one unit per gene spanning its weight window, so GWAS and LD are read per gene) next to the existing region grouping. New optional CV prefilter (--min-cv-rsq,--max-cv-pval) that skips genes where no context/method passes, so their GWAS input is never built.--min-cv-rsq > 0. The CV metric names are configurable (--rsq-option,--rsq-pval-option). Methods without CV metrics (e.g.ensemble) are skipped. Genes with no CV results are kept; genes whose weight file cannot be read are kept with a message.--unit genestops if a gene has more than one xQTL meta row, and makes the unit ID file-name safe.--unit genestops early when a gene has no start/end in the xQTL meta, and when two gene IDs would give the same file name.twas.R: a missing GwasSumStats file still means an empty window, but a missing folder now stops as a wrong path. Thetwasstep removes an old GwasSumStats file before rebuilding it, and makes the gene part of output names file-name safe.--twas_unit,--twas_prefilter,--rsq_optionand--rsq_pval_option.No default changes for existing users:
twas_ctwas.ipynbkeeps the region unit and the prefilter off. Opt in with--twas_unit gene(outputs are then named<name>.<gene>.<gene>.twas.rds) and--twas_prefilter. Use the prefilter for TWAS-only runs, since the cTWAS steps read the same manifest and cTWAS should use every gene.Tested on SCC with pecotmr 0.9.6: MSBB eQTL TWAS on all 22 chromosomes gives the same per-chromosome gene counts as the 0.9.5 run, and the TWAS/cTWAS integration and script tests pass (27 tests).